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I .  INTRODUCTION 


A  previous  report  [3]  presented  a  general  outline  of  procedures 
for  simultaneously  incorporating  various  sources  of  auxiliary  informa¬ 
tion  into  an  impact  acceleration  injury  prediction  model.  This  report, 
which  serves  as  a  companion  volume,  discusses  computational  aspects  of 
the  procedures.  Specifically,  the  report  presents  computational  proce¬ 
dures  which  allow  the  application  of  commonly  used  nonlinear  estimation 
programs  found  in  statistical  packages  such  as  BMDP  [1]  and  SAS  [2]. 

Throughout  this  report,  reference  will  be  made  to  equations  in  the 
previous  report  [3].  These  equations,  identified  as  (1)  through  (17) 
in  that  report,  will  be  referenced  herein  by  the  same  numbers  as  a  matter 
of  convenience.  Equations  in  this  report  will  therefore  begin  with 
reference  number  (18).  In  addition,  the  notation  used  in  the  previous 
report  will  be  adopted  without  repeating  Its  definition.  Therefore,  the 
reader  may  find  it  convenient  (or  perhaps  necessary)  to  have  that  report 
readily  available. 
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II.  COMPUTATIONAL  PROCEDURES 


To  estimate  the  model  In  (17)  by  using  a  statistical  package  such 
as  SAS  or  BMDP,  one  should  iteratively  compute,  for  each  observation: 

(i)  the  model  derivatives  of  (17)  with  respect  to  the 
parameters  and  p 

and  (ii)  the  regression  weights,  [p^(l  -  p^)].* 

The  derivatives  of  (17)  with  respect  to  JL  and  p  are 

2  — ^ 

_i_  Pi  =  4>(U.  (0.,  P))(l  -  p  > 

94 

|_Pi  -  ♦(0i^61,  p)>  ♦  {-(1  -  p2)~h  S±  +  [P/(1  -  p2)]u±(^1,  P)} 

where  4>  (U)  =  d_  0(U). 

dU 

This  paragraph  concludes  the  discussion  of  the  problem  without  regard  to 
prior  information. 

A.  PRIOR  INFORMATION 

The  model  parameters  in  (6)  or  (7),  as  the  case  may  be,  can  be  esti¬ 
mated  conveniently  and  efficiently  by  using  the  model  in  (17)  to  incorporate 
the  prior  information  into  the  parameter  estimates.  Here  again  the  intuition 
behind  this  approach  comes  from  the  likelihood  of  the  observations.  Consider 
the  likelihood  in  (12).  The  form  of  this  likelihood  suggests  the  models: 


(18a) 


yu  -  *[<1  -  p2)"*5^  -  p8i)]  +  C1 
y2i  ‘  ^2  +  £2i 


(18b) 


where  E^)  -  0,  Var  (S^)  =  p1(l  -  p±),  pt  -  ♦  [Uj.flJp  p)] 

E(e21)  =  0,  Var  (e21)  *  o2  . 

Note  that  the  and  e error  terms  are  stochastically  Independent, 
since  yli/si  and  y  are  independent.  This  fact  makes  the  analysis  by 
statistical  packages  convenient  since  their  nonlinear  regression 
algorithms  cannot  analyze  correlated  observations  cirectly. 

As  stated  in  the  previous  report  [3],  it  is  assumed  that  the  prior 
information  can  be  in  the  form  of 


(i)  a  priori  estimates  of  some  function  of  the  model  parameters 

and/or  (ii)  a  priori  knowledge  in  the  form  of  model  parameter  equality 
constraints. 


The  technique  used  to  model  the  prior  information  discussed  in  the 
previous  report  can  be  applied  to  the  models  in  (18a)  and  (18b)  to  form 
an  overall  model  of  direct  and  auxiliary  information. 

Let  the  prior  estimates  of  some  function  of  the  model  parameters 
be  expressed  by  (6),  where 


Recall  that  it  is  assumed  that  the  prior  information  embodied  in  _r  is 
stochastically  independent  of 
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Then  all  the  information,  barring  parameter  constraints,  can  be  modeled 
as 

r  -  &<8,  p)  +  v  (19a) 

L  m  Id*  P)  +  n  (19b) 


where  E(\0  =  0,  Var(\j)  =  ^ 


-  *(u1(81,  p))  . 
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The  vector  r^  is  conceptually  and  computationally  considered  as  an 
observation.  Due  to  the  structure  of  most  statistical  packages,  this 
concept  can  be  computationally  exploited  only  if  ^  is  a  diagonal  matrix, 
where  the  diagonal  elements  of  ^  are  considered  as  weights  for  the  observa¬ 
tions  in  Even  if  jP  is  a  diagonal  matrix,  it  is  probably  the  easiest 
for  the  user  if  he  or  she  transforms  the  elements  of  r_  so  that  the  transformed 
elements  have  the  identity  matrix  as  covariance  matrix.  This  is  done  in 
this  report. 

The  elements  of  r^  can  be  transformed  by  premultiplying  r  by  a  matrix 
T  such  that 

1  _i 

II-I  * 

If  r*  denotes  the  transformed  elements  of  then  (19a)  can  be  reexpressed  as 

r*  *  r_r_  =  £*(§,  p)  +  Tv 
where  !*(§»  P)  “  £ g(3,  p) 

and  E(r v)  *  0,  Var (T v)  =  £  (£  £)_1£*  “  I. 


B.  PARAMETER  CONSTRAINTS 

To  include  constraints  on  the  parameters,  let 

<B,  p)  "  h(e)  (20) 


where 


dim  9_  <_  dim(j5,  p). 


Jt  ■  •'! 


Then,  it  follows  that,  defining  v  -  Tv, 

r*  “  (h  (0)  )  +  V* 
£  =  £(h(0))  +  J1 


(21a) 

(21b) 


Note  that  indicator  variables  can  be  used  to  facilitate  the  computation 
of  the  estimates  of  (j5,  p)  with  a  statistical  package. 

For  example,  consider  the  case  without  constraints  and  let 

zi  “  5liyli  +  ^2iy2i  +  ^  ^  "  52i^r3 

k^g,  P)  »  <Slif1±(l,  P)  +  $2if2i(£*  P>  +  (1  -  6li)(1  -  62 i>8j(!»  P> 


ei  "  6n5i  +  62iEli  +  (1  <Sli)  (1  "  62i)vj 


where  6, 


62i 


i  if  zi  «  y11,  k^g,  p)  =  fli(6,  P),  and  e.  = 
0  otherwise 


(l  if  zf  =  y2i»  ki(l»  p)  *  f2l(§*  and  ei  = 


i 0  otherwise 

and  i  s  1 , . • . , n  ,  j  *  l,...,m  ,  m  d im  r ^  . 

Here  r ^ ,  g*,  and  Vj  are  the  elements  of  r_,  £*,  and  v>  respectively. 
Now  rewrite  (22)  as 

z  -  k(g,  P)  +  e  (23) 


where  z 


,  k(g,  p) 


kx  (jB,  P) 


kn(£,  P) 


III.  A  POSSIBLE  EXPERIMENT 


To  elucidate  the  application  of  estimating  the  model  in  (23),  an 
example  of  a  possible  impact  acceleration  injury  experiment  with  n 
observations  may  be  considered.  Suppose  that  for  i  ■  l,...,n: 


(i) 


x  =  (1,  x,  .)'  where  x  denotes  peak  sled  acceleration  for 
— i  li  li 

the  i^  subject, 


(ii) 


y  is  a  dichotomous  injury  observation  for  i1-*1  subject,  with 
value  of  1  for  injury,  0  for  no  injury. 


th 

and  (iii)  y  is  a  continuous  preinjury  measurement  for  the  i  subject. 
Then  the  model  of  the  empirical  data  as  given  in  (19b)  has 


3  - 


9 


h 


B 

B 


01 

11 


,  B 


2 


3 

3 


02 

12 


A.  SOME  ASSUMPTIONS 


Suppose  further,  that  for  a  particular  level  of  peak  sled  acceleration, 
Xq  say,  a  good  prior  unbiased  estimate  of  injury  probability,  p^,  is  avail¬ 
able.  Let  this  prior  injury  probability  and  peak  sled  acceleration  level 
be  related  by 


*^®01  +  ^11*0^ ’ 


Assume  that  the  variance  of  p^  is  known  or  estimated  to  be  ^  . 

Thus,  the  prior  information  as  modeled  in  (19a)  can  be  written  as 
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r*  =  &*(£,  p)  +  Tv 

where  £  is  a  lxl  matrix  estimated  by  <|>  \  r*  is  a  lxl  vector  equal  to 
'T1P0,  and  £*<£,  p)  =  i|/_1[$([l,  xQ,  0,  0]B)]. 

Finally,  assume  that  the  elements  of  £  are  constrained  by  the  equation 
3^2  =  T**is  constraint  can  arise  quite  naturally  from  the  following 

situation.  Suppose  that  the  event  of  a  preinjury  measurement,  (e.g.,  change 
in  evoked  potential  response)  exceeding  some  critical,  prespecified  value, 
y^,  is  always  at  least  as  likely  as  the  occurrence  of  an  injury,  for  all 
levels  of  peak  sled  acceleration.  This  can  be  mathematically  represented 
by 

pr(y2i  >  yQ|x)  £  Pr^ii  =  M*)  for  a11  21* 

Now, 

Pr<y2i  >  y0!x)  *  $[-(y0  “ 


and  Pr(yli  =  l|x)  =  . 

Therefore,  for  all  x* 


*l-(y0  -  1  1 

However,  this  can  only  be  true  if 

612/02  '  611  • 

The  constraint  imposed  by  the  vector  valued  function  h  can  be  described 


by 


8oi 

r  1 

0 

0 

0 

611 

=  h(6)  - 

0 

1 

0 

0 

802 

0 

0 

1 

0 

S12 

0 

2 

a 

0 

0 

.  0 

0 

0 

0 

1 

0  , 
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where 


However,  depends  upon  a 2,  which  is  unknown.  In  other  words,  the  prior 
information,  in  terms  of  the  parameter  constraints,  is  not  complete.  In 

A 

this  particular  case,  the  best  that  can  be  done  is  to  estimate  h  by  h. 


and  5 2  is  the  minimum  variance, 

preinjury  data,  i.e.  o„  =  — 2 _ 

‘  n  -  2 


unbiased  estimator  of  based  on  the 

A  A 

02,  a2  being  the  maximum  likelihood 


estimate  of  c^. 


B.  MODEL  REPRESENTATION 

The  model  representation  for  computer  analysis  is  then  given  (for 
i  *  1 . n  and  j  =*  1)  by: 

zi  "  6liyli  +  62iy2i  +  (1  -  6li>(1  -  62i)rj 
ki<6,  P>  =  6lifii(E(£)  +  62if2i(^)  +  (1  -  6li)(1  -  52i)8j^(i)) 
ei  “  6ii^i  +  ^2ie2i  +  (1  -  <5u)(l  -  <S21)V* 
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where  6  and  were  defined  previously,  and 


f  !i  (— <— )  > 


P)) 


f2i<M£)) 


^02  +  XilS20ll 


g*(£(0))  -  <rVfi01  +  *0Bu>. 


If  a  separate  estimate  of  B  based  on  the  dichotomous  injury  data 
and  separate  estimates  of  3^  and  0 2  based  on  the  continuous  preinjury 
data  are  available,  then  a  "quick-and-dirty"  estimate  of  p  can  be  computed. 
Note  that 

ECE2iyii^  "  -P<J21M5i§i)  (24) 

Var(e21y1;L)  *  -  p2[  +  ^(xjj^)  ] ) 

The  results  in  (24)  are  derived  in  the  appendix  of  this  and  the  previous 
technical  report  [3].  From  (24)  it  is  easy  to  see  that 


E*~e2iyli<a2<K-i-l))~1*  “  P 


Thus, 


,f1£2iyii/[*<^l)"°2J 


is  an  approximately  unbiased  estimate  of  p.  Here  is  the  probit  estimate 

o 

of  3^  based  on  the  dichotomous  injury  data  only  and  a2  is  the  unbiased 


estimate  of  02, 


n  2 

1  I  . 

(n  -  2)  i-1  21 


Unfortunately,  the  estimate  In  (25)  can  sometimes  have  absolute 
values  exceeding  one.  A  more  refined  (approximately  unbiased)  estimate 


of  p  can  be  easily  obtained  by  using  an  iterative  least  squares  statistical 
algorithm.  To  compute  an  iteratively  reweighted,  least  squares  estimate  of 
p,  simply  let 


e2lyil 


e2nyln 


(26) 


be  the  dependent  observations  and  let 


•32»(x^i) 


(27) 


be  the  independent  observations.  Compute  the  initial  weights  from  the 
reciprocals  of  the  variances  of  C2iyn  picking  an  initial  value  for 
p.  The  model  derivatives  are  estimated  by  (27). 


C.  CORRELATION 

Thus  far ,  all  models  have  assumed  that  p  is  constant  for  all  x  . 

An  approximate  chi-square  test  can  be  conducted  to  test  the  null  hypothesis 
that  p  does  not  depend  upon  x  .  From  (24)  the  means  and  variances  of 
£2iyli  can  computed.  Let 


P1  ’  E<C2lV  ’ 


a  =  Var(e„  y  )  , 

i  2r  li 


e2iy2±  “  pi 


and  w  *  i  I  w 
n  i=l  1 

Then  w  , ...,w  are  stochastically  independent,  each  with  mean  zero  and 
1  n 

variance  equal  to  1.  It  follows  that 

n  -  2 
I  (w  -  w)Z 

i-1  1 

has  an  approximate  chi-square  distribution  with  (n  -  1)  degrees  of  freedom. 

Of  course  all  the  y^'s  an^  ai*s  are  not  known»  but  they  can  be  estimated  by 
substituting  the  estimates  of  the  model  parameters  via  the  equations  in  (24). 
Replacing  the  y^'s  and  cr^'s  by  their  estimates  should  still  yield  a  reason¬ 
ably  powerful  test  against  significant  departures  from  the  null  hypothesis 
for  moderate  to  large  samples  sizes. 

If  the  chi-square  test  rejects  the  null  hypothesis,  then  it  is  still 
possible  to  estimate  a  probit  model  where  p  depends  upon  x  .  However  in  this 
case,  it  would  be  desirable  to  increase  the  sample  size.  This  is  due  to 
the  addition  of  the  extra  parameters  used  to  model  p  .  The  computation  of 
the  model  parameter  estimates  would  also  be  considerably  more  difficult  to 
obtain. 
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IV.  DISCUSSION 


The  purpose  of  using  auxiliary  information  is  to  improve  the  estimate 
of  3^.  In  some  estimation  situations,  however,  the  use  of  auxiliary  infor¬ 
mation  may  result  in  estimates  of  j5  that  are  slightly  less  accurate  than 
if  the  auxiliary  information  were  not  used  at  all.  If  the  constraint 
information  is  correct  and  the  prior  estimates  of  3^  are  close  to  the  true 
value  of  3^,  then  the  a  priori  auxiliary  information  should  always  contribute 
to  reducing  the  mean  square  error  (MSE)  of  _3  . 

If,  however,  the  sample  size  is  small,  there  is  no  apriori  Information 
linking  the  parameters  of  3^  with  3^»  and  in  addition  the  correlation  between 
the  injury  tolerance  and  the  side  effect  (i.e.,  preinjury)  is  low,  then  the 
incorporation  of  preinjury  data  could  possibly  contribute  to  a  slight  increase 
in  the  MSE  of  the  estimate  of  3,^.  In  order  to  assess  the  benefits  of  the 
inclusion  of  preinjury  data  for  small  sample  sizes,  a  Monte  Carlo  study  is 
being  conducted.  This  study  will  be  discussed  in  a  future  report. 
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APPENDIX 


_ rr\i 


Since 

2  2 
Var(e21y1;l)  -  E{<E2iyli^  >  -  » 


it  is  sufficient,  in  view  of  the  appendix  in  [3],  to  derive  an  expression 
for  E{(e2ly11)2}. 

Now, 

E{(e21yii)2}  “  £tE2iyli^ 

00  00 

=  ^  e2iIi^tiJfe.T<e2i*  tl)de2idti 


■f0  '  E2i^(E2i|tl>£T<tl>dEHdtl 

.00  .00 


■  /„EU2i|ti)fT(ti>dtl 


However , 

E{£2i|l:i}  “  Var(e2ilti)  +  {E<e2i|,:i)}2 

2  2  2  2  2  2 

-  (1  -  pz)o2  +  pzaz(ti  -  pT  )VaT 

since  e2ilti  has  a  normal  distribution  with  mean 

pa2(ti  -  v  )/aT 


and  variance 

(1  -  p2)02  . 

Thus, 


-16- 


A 


E^e2iyu)  }  -  (1  -  P2)o2/  fT<ti)dt1  +  (P2a2/o2)/_ao(t1  -  yT  KMt^dt 


Now,  from  (3)  and  (4)  in  the  appendix  of  [3], 
0  “UT  /0T  1  -u2/2 


*(-p  /a  ) 
Ti  T 


we,)  . 

~i~ 1 


Also, 


-UT  /°T 

0  2  2  i  :  1  2/0 

/  t  Ct  -  y  )  la41] f  (t  )dt  -  /  u2  /5?  e“U  /2  du  . 

_oo  1  Xi  T  Ai  1  i  _00 


By  performing  integration  by  parts,  the  integral  on  the  right  hand  side 
becomes 


-u(2tt)  e 


>3  -V 


-y  /a  -y  /ct_ 

1  1  -i.  ->iu2 

+  /  7^  e  du 


This  expression  is  equal  to 

-<yT  la  )0(-y  la  )  +  4>(y  la  )  -  J^CxlS. )  +  <Kx!31) 

I  ^  X  Ti  T  — 1 — A  — 1 — 1  — 1—1 

Therefore, 

E^e2iyli)2}  "  (1  "  p2)a2l,,^-i^i)  + 


and  hence, 


Var(e2iyli)  "  a2[*(^l) 


p2((xie1)<t(x^B1)  +  02(x|61))] 
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